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The Szekeres inhomogeneous models can be used to model the true lumpy universe that we 

observe. This family of exact solutions to Einstein's equations was originally derived with a general 

metric that has no symmetries. In this work, we develop and use a framework to integrate the 

angular diameter and luminosity distances in the general Szekeres models. We use the afhne null 

geodesic equations in order to derive a set of first-order ordinary differential equations that can 

be integrated numerically to calculate the partial derivatives of the null vector components. These 

equations allow the integration in all generality of the distances in the Szekeres models and some 

c7^ . examples are given. The redshift is determined from simultaneous integration of the null geodesic 

Cn ' equations. This work does not assume spherical or axial symmetry, and the results will be useful 

►;^._^' for comparisons of the general Szekeres inhomogeneous models to current and future cosmological 

jrt ■ data. 

PACS numbers: 98.80.Es,95.36.+x,98.80.-k 

I. INTRODUCTION 

\m^ ' Recently, there has been a renewed interest in studying some current questions in cosmology using inhomogeneous 

\^ . models, see for example the reviews [l|, Q and references therein. Studies have included for example the formation of 

structures in the universe [j-Q and the question of cosmic acceleration, e.g. see [7Hl8{ and references there. It is also 

Q,. of interest that such inhomogeneous models can offer a wider range of interpretation of cosmological observations. 

I ' Whereas one can find a large body of literature studying exact solutions to Einstein's field equations that are 

2 [ inhomogeneous cosmological models [20, [2l| , relatively little work has been done comparing inhomogeneous models 

to various observations. 

jrt ' In this paper, we consider the calculation of the area and luminosity distances in the general Szekeres models. The 

affinely parameterized null geodesic equations are used to derive a set of first-order ordinary differential equations 
ready to be integrated numerically to calculate the partial derivatives of the null vector components. The results allow 
the integration of the distances in the general Szekeres models. The redshift is calculated from numerical integration 
Qs of the set of null geodesies. The general results should be of immediate application to comparisons between the general 

QQ ] Szekeres models and cosmological data. 

Q^ ' The models were originally derived by Szekeres [22, [2j] as an exact solution to Einstein's equations with a general 

OQ ' metric that has no symmetries (Le. no Killing vector fields) and a dust source. The models have been investigated 
If^ , analytically by several authors [S, [a, [25l - [3ll [3ll - l34l |37| . Reference (35| classifies them in the same category as the 
^^ ■ observed lumpy universe (see the table on page 37 there) and the models are considered as one of the best exact 
^^ i solution candidates to represent the true lumpy universe we live in. The models have been discussed in regards to 
^~^ •■ structure formations in [a, |6|, |28| and also most recently in the context of cosmic acceleration and comparisons to 
J>. [ supernova data in [l6|, [l9j and the Rees-Sciama effect in the cosmic microwave background radiation in ^] . 
• '~j ■ It is perhaps worth clarifying that we are not proposing the Szekeres model as the true and ultimate model of the 

K% ' universe but rather developing a framework based on exact inhomogeneous cosmological models where cosmological 
S . observations can be investigated with a wider range of possible interpretations. 

II. THE SZEKERES MODELS 

The models have been written in at least three different sets of coordinates: The original set was given by Szekeres 
in [22, 23], then Goode and Wainwright [2^ proposed a second set of coordinates, and then finally a third set of 
coordinates was used by [a, la, 1221; which we adopt in this paper. The models have been discussed in great detail in 
[20I I24I I33I l34| and we give here only a brief introduction to set the notation. 
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The Szekeres metric can be written as 



ds^ = -di^ 



e — k{r) 



R^ 

E2' 



dr^ + tl^{dp^ + dq^) 



(1) 



where R — R{t, r) is the areal "radius". The function k{r) is related to the energy per unit mass and determines the 
curvature of the spatial sections t = constant. This function divides the models into sub-cases: hyperbolic (fc(r) < 0), 
parabolic (fc(r) = 0), and elliptic (fc(r) > 0). The function E ~ E{r,p,q) is given by 



E{r,p,q) 



(p2 + q2) P{r) Q{t) 



2S{r) S{rY 

and the functions P(r), S{r), Q{r), and C{r) satisfy the relation 

PHr) Q^r) 



C{r) 



S{r) 



Sir) 



C'ir) 



2S{r) 2S{r) 



(2) 



(3) 



but are otherwise arbitrary. 

The geometrical constant e determines whether the {p,q) 2-surfaces are spherical (e — +1), pseudo-spherical (e = 
— 1) or planar (e = 0). That is, the constant e determines how the 2-surfaces of constant r foliate the 3-dimensional 
spatial sections of constant t. The function E{r^p, q) determines how the coordinates (p, q) are mapped onto the unit 
2-sphere, pseudo-sphere or plane for each value of r, see for example 24, 32, 33]. 

The Einstein field equations with a dust source and no cosmological constant read 



{Rtit,r)f 



2M{r) 
R{t,r) 



k{r) 



and 



8Trp{t,r,p,q) = 



2(M,,(r)-3Af(r)%:) 
R'^{R,r-R^) 



(4) 



(5) 



where we have set G = c = 1 and the function M{r) represents the total active gravitational mass in the case e = +1 
[32,l33]- The evolution of R{t, r) depends on k{r) and is given as follows: 
The hyperbolic case: k{r) < 



R(t,r) 



M{r) 
-k{r) 



(cosh?7 — 1) 



(6) 



The parabolic case: k{r) — 



The elliptic case: fc(r) > 



t- tsir) = a 



M{r) 



(-fc(r))3/2 



(sinh r] — rj) 



R{t,r) =M{r) 



^- 



t~tB{r) = aM{r)^ 
6 



i?(i,r) = ^1^(1 -cost;) 



k{r) 



(7) 

(8) 
(9) 

(10) 



t-tB{r) 



cr --» (^ - sm rj) 



{k{r)f'^ 



(11) 



where tsir) is an arbitrary function of r and represents the Big Bang or Crunch time. Also, u = ±1 is in order 
to allow for time reversal. There are several sub-cases of the Szekeres models, and a given model is specified by six 
functions that can be reduced to five by using the coordinate freedom in r. This indicates their rich geometry |32h34|. 



III. A CONVENIENT FORM FOR THE NULL GEODESIC EQUATIONS 

The null geodesic equations govern the propagation of light rays in a given spacetime and are necessary to solve 
in order to derive observable functions for a given model. The non-affinely parameterized null geodesic equations are 
given by 

where r is an arbitrary (non-affine) parameter, A(r) is a function of r, and k°' = ^^ is a null tangent vector to the 
geodesies. Now, let's recall that if one applies the following change of parameter (see for example [23]) for a discussion) 



T -^ s{t) = / Yl--dT 
Jro Hz) 



(13) 



where C is a constant and it follows in this parameterization that 

k°';fik^=0 (14) 

where the null tangent vector fc" — ^- and the geodesic equations ([T4l) are all affinely parameterized. 

The afhnely parameterized null geodesic equations for the Szekeres model are given by [2, ST], (see footnote |38|): 



' R,r-R^ V i?,r-i?% 2(e-fc) 



R E,rE,p EE,^, ^^^^ ^ ^R E,rE,g EE^^^ ^^^^ _ R e k , ^ ^ ^ ^ 

£^2 R,r-R^ E^ R,r-R^ E^R,r-R^ 

kv + 2:|ifc*F - l:^^!!_^(ij,, ij,p -ii;i?„. )(fc'-)2 + 2 (^^ - :|!:^ ff 
_:^(P)2_2:^Ffc9 + :^(fc9)2 = o (17) 

i? E E 

ki + 2^fc*fc^ - ^ ^';"^^"^ (i^,. i?„ -EE,,r wr + 2 (^ - %^) fc'-fc' 

+:^(fcp)2_2:^fcpfc9_:^(fc9)2 = (18) 

E E E 

where ' = ^ and we will use the two notations interchangeably. Now, in order to rewrite these equations in a 
convenient form for us to solve for the null vector components, we proceed in this way. First, we define the compact 
functions 

H = ^'^''- ' ^^f^^' (19) 



and 



The geodesic equation (|15p then becomes 



F^^. (20) 



fc* + \{H)Akl' + \{F^)A{k''? + [k'^f] = 0. (21) 

We multiply equation (|16p by H and use the total derivative expression 

^ = ^^^k^ + ^ ,.r + *,pfcP + 1',,/c9 (22) 



in order to rewrite the second geodesic equation p6|) as 



Hk- + Hk- - \{H)Ak''? ~ \{F^),r[ik''? + [k^f] = 0. 



(23) 



2' '■' ' 2 
Similarly, we multiply Eqs. p7|) and (|18p by F^ and use ((22)) to rewrite the third and fourth geodesic equations as 



F^F - i(iJ),p(fc'-)' + {F^)k^ " ^(i^'),p[(fc'')' + [k'^f] = 



(24) 



i^'fc^ - IwAk^-f + (^')- fc« - ^(i^').,[(fc^)' + (fc')'] = 0. 

A relation that will be useful for simplifications is the null vector condition k°'ka = and reads 



(fc*)2 - H{¥f - F^[{kPf + {k^f 



0. 



(25) 



(26) 



Now, the null geodesic equations (|2ip . ([23]), (12^ . and ([25]) constitute a system of 4 second-order ordinary differential 
equations (ODEs) for the functions {i(s), r(s),p(s),(7(s)} where the coefficients are composed of the metric functions 
evaluated on the null cone using the field equations and the model specifications. We use a fourth-order Runge-Kutta 
algorithm with adaptive step size for the numerical integration [33| . The code iterates between calls to evaluate the field 
equations on the null cone and calls to integrate the ODEs. We implement the Runge-Kutta code with the function 

vectors (see for this notation M) given by y = {i,r,p,g, ^, £, % g} and ^ = {^, £, g, g, 0, 0, |f , ^} so 
our system of 4 second-order ODEs is transformed into a system of 8 first-order ODEs. 

While this integration provides us with the four components {fc*, fc'', fc^, fc''} and is enough to compute the redshift, 
we need to solve for the partial derivatives of these components in order to calculate the area and luminosity distances 
and we do that in the next two sections. 



IV. A SET OF EQUATIONS FOR PARTIAL DERIVATIVES OF THE NULL VECTOR COMPONENTS 

First, we use the null vector condition (^5)) into the r-component of the null geodesic equation (p5| to write 



HJ-{¥) + ¥^^H-\{H)Ak 



'''' ~l^-^[ikr-Hikr]=o. 



Next, we note that we can use (|22|) to write the useful relations 
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as 



= j^{{F'),,k^)^iF'),pak^ + {F'),pk%. 



dx' 



Now, taking the partial derivative of equation (P7|) wrt x" and using equations (P5| and ([2^ . we get 



H 



I'iKa) + Kpk „ 



ds 



H,c^{kn + [H,0c.k^ + H^pk%] F + H^pk^V^ - ]-H^ra{K'f - H^^k^V, 



ds 



1 f{F^). 



2\ F^ 



.) {{kr H{¥r) - i^ [k'k]^ - \HAkl' - Hk^k^^j = 



(29) 



(30) 



(31) 



where -^{k^) is given from equation ([23 



i^'-^-'n 



k''i:H-\{H)A^r - \^^W? - H{k^f] 



d 

ds 



(32) 



In a similar way, we take the partial derivative of the i-component of the null geodesic equations, (|2ip. and use the 
relation ([25]) to obtain 



d_ 

ds 



ik%) + klpk% + ^HMkr + Htk^k:, + i (^^) {{kr - H{¥f) 



(P'),^l,j.t 



F2 



k'k' 



1 



H,4k^Y-Hk^k^^] -0. 



(33) 



Similarly, from the partial derivatives of the p-component and (/-component of null geodesic equations, i.e. (1241) and 
d^S]) . and using the relation ([50]) . we find 



ds 



y^^a) + ^.13^. a 



+ (^'),a^(fc') - ^^,P«(fc')' - H,,k-¥^ + [{F^),t3o.k^ + {F^),pk%] F + 



{{F'),f>k^) k^. - 1 (^) {{kr - H{kr) 



^-^Uki^-\HAki'-Hk^k:^ 



and 



ds 



\kn) + k ak 



+ {F^),a-^ik') - Ji?,ga(fc^)' - i^.^fc^A:;-, + [iF'),0^k^ + iF'),pk%] fc« + 



ds 



{{F^),,kP) ki - 1 ^^^ {{ky - ff(fc'^)2) 



(^'),. 



p2 y- %a 2 
Finally, we apply the same process to the null vector condition (j26p to write. 



H^^{k-)^-Hk-k-^] =0 



(34) 



(35) 



fc': = 



1 

2"F 



-^(kn' + ^ {2k'k% - (F2).„ {{kPf + (fc^)2) _ 2F2 (fcPfcP, + fc^fc^^J) 



(36) 



Now, equations (PT|) . (1551) . ([M| . and ([55]) provide a set of 16 first-order ordinary differential equations (4 equations 
for each a — t, r,p, q) that can be integrated numerically for the 16 components {fc*^,, k^^, k^^, fc^^}. This integration is 
done simultaneously with that of the system from the previous section for the 8 first-order ODEs for the components 
fc*, k"^, k^, and k'^. The system of the 16 ODEs here is given in the appendix as equations (A1)-(A16). 

From a numerical point of view, we found it more practical to solve the system of the 12 ODEs given by equations 
(I55)) . ([Ml) . ([55)) plus the 4 equations given by ([55| (that is the set of equations (A5)-(A20) given in the appendix). 
Again, we integrate the system of ODEs using a fourth-order Runge-Kutta algorithm with adaptive step size |39|. 
The code iterates between calls to evaluate the field equations on the null cone and calls to integrate the ODEs. 

With the integration of the null vector components as well as their partial derivatives, we can now calculate the 
area distance, the redshift, and the luminosity distance. 



V. THE AREA AND LUMINOSITY DISTANCES FOR THE SZEKERES MODELS 



As usual, the area distance. Da, is related to the surface area, SS, of a propagating light front of a bundle of light 
rays by the relation (see for example f24|) 



SS = D\6n 



(37) 



where 5il is a solid angle element. 

Now using the relation of the area distance to the expansion optical scalar 9 (see early work by [40, |41| and also 
[2J| for a recent discussion), one obtains 



d\n{SS)^20ds 



(38) 



where s is an afSne parameter and 6 is given in terms of the afBnely parameterized tangent vector as 



-2^";"- 



It foUows from the three above equations that 
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The right-hand side integrand of equation (PH)) can be evaluated as 
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where the connection coefficients can be obtained straightforwardly from the metric as: 
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Putting these into equation PO]) gives 
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(42) 



and upon applying the total derivative equation (P^ . we get the following equations for the area distance: 



Da = F H^/'^ exp 



(fc;t + fc;; + Fp + fc«jds 



The luminosity distance is given from its relation to the area distance expression (j43p as 



(43) 



(44) 



Given explicit functions for a Szekeres model, the area distance is calculated from equation (|43p where the components 
fc*(, fc''^, fc^, and k'^„ are numerically integrated from the system of ODEs derived in section IV and listed in the 
Appendix. The redshift is also numerically integrated from the null geodesic equations as described in section VI 
further. 
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FIG. 1: Luminosity distances for a Szekeres model that is not axially or spherically symmetric. To the left, the value of g is 
fixed to —200 while p is varied by taking the values —100, —50, 0, 50, 100. To the right, the value of p is fixed to —100 while q is 
varied by taking the values —200, —100, 0, 100, 200. The Szekeres inhomogeneous model used here is for an illustration purpose 
only and was introduced in ^5]] within structure formation context. The model is specified in our section V-A. The luminosity 
distance for an open FLRW model is plotted as well. 



A. Examples of luminosity distance versus redshift plots 

For illustration, we integrate and plot luminosity distances versus redshift for a Szekeres model that is not axially 
or spherically symmetric. We use a model that was introduced in [5|] within large scale structure formation context 
and we simply use it here for illustration purposes. In this case, we set: 

. e = +l. 

• is iy) — for a simultaneous Big Bang. 

• M(r) = (sinh(r))3. 



k{r) 



(i.e. at large r, the spatial curvature goes to zero as indicated by CMB observations [42|) and at 



very small r curvature goes negative in accord with local observations of matter abundances with lower density 
than the critical density. 

• The functions {S, P, Q} are given by model-1 of ^ with {140, 10, -113 ln{l + r)}. 

The results for various values of p and q values are given in Figure I along with an open FLRW model. The exploration 
of more sophisticated Szekeres models and comparisons to supernova data and other cosmological distances data will 
be presented in follow-up investigations. 

B. Application and verification of special cases 



We derived a general expression for the area distance and luminosity distance for the Szekeres models, 
important to verify that these expressions reduce to ones that we know in the special cases. 



It is 



1. Axially symmetric case 
As explored in [2|, |37| , the conditions on the Szekeres metric functions for axial symmetry are given by 






(45) 

EE^qr = E^qEr. (46) 

It is trivial to see from our compact notation in the geodesic equations (|24p and (1251) that these axial symmetry 
conditions are here simply 

if,p = H^q = 0. (47) 

Conditions (H71) assure that dq = dp = holds along the whole geodesic. We proceed then by setting dp = dq = 
into equation (pG)) so (fc*)^ — H{k^Y = or 



t\2 



(kr = j^{k') 



Using this in equation ((2T|) and integrating for (fc*)^ gives 



(fcT=exp(-/^*^4 



Now, from (HSl) . we have 






H 



and using it in equation (j49]) and taking its square root gives 

/c* = exp ( / {VH),tdr\ 



(48) 



(49) 



(50) 



(51) 



1 ff,i 



where we have also noted that 2~/h ~ (v ^),t- 
Next, using (P5|) . we get immediately for fc'': 



k'=-^cxpfl\VH),tdr 



These two expressions agree with the ones obtained in 2] for the axially symmetric case. 
Now, our expression (j42|) for the area distance reduces to 



dlnDA = 



2 V F2 F2 






H 



t T" Kr) 



ds 



(52) 



(53) 



where we use the tangent vector component equations (j5ip and (1521) for axial symmetry to express the partial deriva- 
tives in the last two terms in (155)) as 



H, 



A:': 



2^0 H 
2 i? 



-dt 



2 7*0 ^ 



di 



fc''. 



(54) 
(55) 



and where we have used ^dt — — j^dt and ^-7= = {VH)^t in equation ([55|) . We also note that with axial symmetry 
p and q are being constant and we can combine the term in ()54p plus the last term in ()55p to make a full derivative. 
Putting this into the integration of the last part of equation ([SS)) gives 
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1 r ii^r 


1 /■* Hi 




/ (fc*t + A:;;)ds = 

/so 
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4X0 H - 


4/0 ^ 



^dt. 



(56) 



These two terms cancel with the integrated two middle terms in equation (|53p yielding a straightforward integration 
result 



Da^F. 
This result is in agreement with the result of [1^ |4^ of the axially symmetric case. 



(57) 



2. Spherically symmetric case 

In the more special case of spherical symmetry, i.e. e — +1 and E^^ — 0, our expression (|42|) for the area distance 
reduces to simply 

d\uD^ = \ i^-^k' + ^-^kA ds. (58) 



2\ R^ i?2 

Integrating both sides gives 

Da = R (59) 

which is the well-known result for the observer area distance for an observer located at the center of the Lemaitre- 
Tolman-Bondi spherical model, see for example |24| . 

VI. THE REDSHIFT IN THE SZEKERES MODELS 

In order to plot the luminosity distance as a function of the redshift, the latter needs also to be integrated in all 
generality. We start with the standard relation 

[k Uaje [k )e ,„„^ 

where fc" is the affinely parameterized null vector, u" = (1,0,0,0) is the 4-velocity vector, and the subscripts e and 
o are for emitted (at the source) and observed (at the observer) respectively. Now, we use equation (PT|) and write 

^I^ = -iH),,k\kr - i^k^kr - H{¥f] = 2fc*^ (61) 



which we use along with taking the natural log of both sides of (|60p and differentiating both sides wrt to the affine 
parameter s to obtain 



1 dz 1 dk* 1 



1 + z ds k* ds 2 fc 



{H)Akn' + ^-^{ikr-H{kn') 



(62) 



While it is informative to see this expression, which one can integrate simultaneously with the three other null 
geodesic equations (l23t . ((24)) . and (|25p to find the redshift, one can equally use equation ((60)) and get fc* by integrating 
the four null geodesic equations. We proceeded with the latter method for our numerical calculations. An alternative 
way to find the redshift and the related drift effects in the Szekeres models can be found in p3| . 

As an aside, we note that using the null condition k^ka = to substitute (fc*)^ — H{k^Y by ^^[(fcP)^ + (fc^)^], the 
redshift equation (|62p can be written as 

iT^S = -^ mAkl' + iF^iikn' + ik^n (63) 

which is precisely the redshift relation that was derived in our previous work [161 144| . 

A. Application and verification of special cases 

We derived a general expression for the redshift for the Szekeres models. It is important to verify that this expression 
reduces to the ones in the known special cases. 

1. Axially symmetric case 

In their recent book y , the authors derived an expression for the redshift for the axially symmetric Szekeres models 
that reads: 

ln{l + z)= dr ■*'• "jI^ . (64) 
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Taking the derivative of this equation with respect to the afBne parameter, s, we get 



1 d{l + z) _ Rtr-R.t^dr 



1 + z ds a/1 - k ds 

Using our equation for the null condition (pS)) and k^ ^ k^ = Q gives 

k"- _ v/r^=^ 

Now we multiply equation (p5|) by -j/r^ = p^ to write 



E , 
- rir 

(65) 



(66) 



l + z ds VT~k k^ ' 

We then use equation (|66p to substitute for fc'' in the denominator of equation (j67p and expand the numerator to get 



1 d(l + Z) _ 1 R^trR^r + RRA^f - i.R,rR,t + RR,tr)^ ,, ,,2 



1 + z ds A:* 1- k 



^,'■^2 (-73 73 I 00 ^-E.' 

(F)2 (68) 



Now, we go back to our expression for the redshift, i.e. equation (|63l) and we put there /c^ = fc'? = as well as the 
definitions of H and F to re- write it as follows 

1 d(l + Z) _ 1 R,trR,r + RR,t{^? - {R.rR.t + RR,tr)^ ,. ,,3 



1 + z ds fc* 1 - fc 



{ky (69) 



As one can see our expression for the redshift when reduced to the axially symmetric case agrees with the redshift 
given in the work of 2] for this particular case. 

2. Spherically symmetric case 

The spherically symmetric case happens in the more particular case where E^r — 0. Substituting this condition 
into equation (p^ gives 

Using equation (I55|) into equation ([70)) gives 

1 rf(l + z) _ i?,t, ,^ 



l + z ds VT^^ 

that upon integration yields 



(71) 



which is the usual result for the spherically symmetric Lemaitre-Tolman-Bondi models, see for example [24| . 

VII. CONCLUSION 

We derived and used a framework to integrate the area and luminosity distances in the general Szekeres models. We 
used the general affinely parameterized null geodesic equations in order to derive a set of first-order ordinary differential 
equations that can be integrated numerically to calculate the partial derivatives of the null vector components. These 
equations allow the numerical integration of the area and luminosity distances in the general Szekeres models. We 
determined the redshift from simultaneous integration of the null geodesic equations. This work does not assume 
spherical or axial symmetry and will be useful for comparisons of the general Szekeres inhomogeneous models to 
current and future cosmological data. 
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Appendix A: Set of ODEs for the partial derivatives of the null vector components 



Equation (PT|) from section IV expands into the following 4 ODEs: 



H 



~r\^,t) + '^,t^,t + '^,r'^,t + '^,p'^,t + '^,q'^,t 



ds 



as 



H^tk^k''^ + H^rk''k'lt + H^pkPk';^ + H^fc^/c^t - -H^rt{k''f - H^rk^'k'^t 



1 {iF'),r 
2\ F^ 



^ ^{{kr-H(kr) 



(F2 



F2 



k kf 



H^tikT ~H¥¥t =0 



(Al) 



H 



ds 



(^/C^J -p /C^/C^. -p tZ ^rZ ^. ~\- K ^tZ ^ ~\- /C^/C^ 



ds 



H,tk'¥^ + H.r¥¥^ + H,pkP¥^ + H^qkl^^ - ]-H^rr{¥f - H^r¥¥r 



1 f{F'), 

2 [ F2 



{ikr-H{¥f) 



F2 



/C fv . 



-H,r{¥f-H¥¥^. I =0 



(A2) 



H 



~ry^,p) + ^,t^,P ^ ^,rk^p + ''^.p^^p + ''^.g^'^.p 



ds 



d_ 
' ds ' 



i?,px(fc'") 



+ [H^tpk' + i?,rpfc" + H^pp¥ + i/,,pA:9 + i/,tA:*p + H^rk'^p + i?,pFp + H^qk^] ¥ 



-H^p{¥Y - H¥¥p = 



(A3) 



^4^^^ 



_rz — — (/c) 4" k fk -\- tz k, -T rc/c 4" f£ k 

^Q ^ ly' i^- ly -' -y i/^ :y ly ly 

+ [i?,t,fc* + iJ^.qfc'- + H^pqkP + H,qq¥ + H.t^q + H^,¥g + H^pkPq + H,qk%] ¥ + 

H^tk^¥q + H,r¥¥g + H.p¥k\ + H.q¥¥^ - ]-H^rq{¥f - H^r¥k\ 

~\ (^) m' H{¥f) ^ (fc^fc*, - lH,i¥r - ijfc'^fc:,) = o (A4) 

where 4^{¥) is written in terms of metric functions and null vector components from equation (j27p (i.e. p2p). 
Equation (1551) from section IV expands into the following 4 ODEs: 

^(fc*J + fc^.fc^, + ¥^¥, + A;*pfcf; + ¥^¥, + \H,uikr + H,¥¥, + \ {^^) {{k'T - H[¥f) 



-^ (k% - ii/,t(fc'-)2 - H¥¥^ = 



(A5) 



12 



d 
ds 



(fc^,) + k'^.k'^, + k'^X, + A;*pF, + klX, + ^H^trik-'f + i/,tFF, + i (^^) ((fc*)' - H{Vf) 



(F2 



i^2 



'* ' k^k\ - -H^r{k''f - Hk^'K'^^ 



^^{k%) + kl,k% + k],v^ + fcfpFp + fc*^fc;'p + \H^,,{k^f + ij,,Ffc;^^ + i (^) ((fc*)' - i^(fc'-)') 



iF').t /,,,, 1 



F2 



k'k'^^--H^j,{ky-Hk-¥\ =0 



(A6) 



(A7) 






' p2 y ,q 2' 

Equation (IMl) from section IV expands into the following 4 ODEs: 



F 



as 2 ' 



-ll'l.P'; J- hP h'^ J- hP h^ 4- hP h-P -^ hP h'* 

[iF^),uk' + {F^)^rtk'- + {F^),ptkP + (F^)^,,^' + (F^)^,^*, + (F^),,^, + iF^),pk'^ + (F^),^^',] F + 

r\2\ 

,t 



{{F'),k' + {F^),x + (f2),,f + (F2),,fc«) fc;; - 1 (i^^ ^ {{kr H{k^r 



M^(k%-lH,{kr-Hk^k:, 







(A8) 



(A9) 



-T~ yfv J, j H- rZ ^rZ J. -r /C^/C^ -f- K „K ^ H- n, ^n, r^ 



ds 



(f2),,^(F) - ii/,p.(fc'')2 - ii^pFfc;; + 



[(F2),,,fc* + (F2)_^^fcr + (F2),p,fcP + (F2),,,fc9 + (F2)^,fc*, + (F2),,fc;, + (F2)^pfcP + [F^) ^qk^] kP + 



((F2),fc* + {F-)A^ + (F2),,F + (F2),,fc^) P, - 1 (^i^) {{kr 



r\2\ 



H{¥) 



iF').P 



r.9 t k'k'^ - T;Hrik''f - Hk'^k'', ] = 



(AlO) 



F^ 



(hP \ M k-P h^ 4- hP h^ J- k-P hP 4- k-P hi 



ds 



+ (^'),p^(fc') - ^^.pp(fc^)' - ^,pfc''^:p 



U'-\2\ 



[(F2),,pA:* + (F2),,,pF + {F^),ppkP + {F^).,,pk^ + {F^),tk% + {F^)A:p + {F^).pk% + {F^),,kl] F + 
7^2 I .p ~ l^Fl^p{k ) ~ Hk kp 



(All) 



F 



(UP \ 4_ UPl-t thPh.^ I UP hP 4- k-P hi 



ds 



+ (^ ),.^(fc') - ^H^p^iky - Hpk^k:^ 

[{F^),tqk' + {F^).rqk- + (F^).pqkP + {F^),qqk'i + {F^) ^tk'^^ + {F^) A^q + (F'),P^,^ + (F'),,^',] ^^ + 

{{F^),k' + (F2),fc'- + (F2),,F + (F2),,fc^) F, - 1 (^i^) ((fc*)2 - iJ(F)2) 
(F'lp f,J,J 1 



F2 



k'k'^-H^,{k^Y-Hk-k\\ =0 



(A12) 
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where j^{k^) is written in terms of the metric functions and the null vector components from equation (p4|) . 
Equation (1551) from section IV expands into the following 4 ODEs: 



F' 



d 



(h.1 ^ Jl_ hi hi 4. hi h"^ 4. J.9 hP _L hi hi 



ds 



+ (F')^,-(k^)--H^,(kn'~H,k^k^, 



[{F')^uk' + {F^).,rtk'- + [F^).ptk^ + {F^).,tk' + {F^).tk\ + {F^),rk:^ + (F^),^^, + {F^),,k'^,] k" + 
{{F'),k^ + {F').rk^ + (F2),,F + (F').,k^) k\ ~ \ ^^^ ^ {(kr - Hik^) 



(F2) 



-^lk%^-H,,ikn'-Hk^k:,]=0 



(A13) 



[{F^),trk' + {F^),rrk^ + (F'),p.F + {F^).^,rk'' + (F^),,^*, + (F^),,^, + (F2),pfcP + (F2)^,fc«J k" + 
{{F^),te + (F^),.^'- + (F2),,fcP + (F2),,fc^) A:^,, - \ f ^"j ((fc*)^ - iJ(fc'r) 



2 V F2 



^ f fc*fc* 



-i?,r(fc 



r\2 



i/rfc" 



(A14) 



F^ 



ds 



(ki) + kik% + k'^.k:^ + k%ki + fc:',fc9. 



+ (^').p^(fc') - ^^..p(fcn' - H^,k-k:^ + 



[(F2),tpfc* + (F2),rpfc'' + {F^),ppkP + {F^)^„k^ + {F^),tk\, + {F')^rKp + {F^),pkP^ + {F^),akl] k" + 
((F^),fc* + (f2),F + (f2),,F + (F2).,fc«) fc^, - 1 (^^) {{kr - H{kr) 



(F'),, 



F2 



k^k'^p - ]-H^p{¥f - Hk'k'-A = 



1 htht 



(A15) 



F" 



— (h1\4. hi ht -^ hi l-r 1 hQ hP I 1,9 uq 



2^ ;.t , /772^ ,r 



[{F'),t,k' + {F^).,rqk^ + {F^).pqkP + (F2),,,fc« + (f2),,A:*^ + (F^),.^, + (F2),,fcP + (F2)^,fc«^] fc^ + 
((F2),,fc* + (F2)^.F + (F2),pF + (F2),,A:'') k% - \ f ^"j ((fc*)^ - H{¥f 



2 V F2 



F2 



1 



^ ' A:*fc*, - -Hjk^ - Hk-k^g ] = 0. 



(A16) 



where ■4-{k'^) is written in terms of the metric functions and the null vector components from equation (j25p . 
Finally, equation p6p from section IV expands into the following 4 equations: 



2fc'- 



^{ky + ^ {2k% - (F^), ((F)2 + (A:'')2) - 2F2 (FF, + fc^A:,',)) 



(A17) 



/v ^ — 



fc'l 



2fc'- 



1 



H, 



'" ('t'')" + ^ (2fc*A;,V - (^'),r {{kn^ + (k^f) - 2F2 (Ffcj; + /fc^fc^^)) 



iJ 



P 2fc'- 



'•« 2fc'- 



_^(F)2 + 1 (2fc*fc;p - (f2),p ((F)2 + (fc^)2) _ 2F2 (FfcP + kVp)) 



-^{kl' + ^ {2k'k]^ - (F^),, ((F)2 + {k'^f) - 2F' {kPkP^ + Ffc^J) 



H 



H 



(A18) 
(A19) 
(A20) 
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